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Abstract. This paper deals with the structural analysis problem of dy- 
namic lumped process high-index DAE models. We consider two methods 
for index reduction of such models by differentiation: Pryce's method and 
the symbolic differential elimination algorithm rifsimp. Discussion and 
comparison of these methods are given via a class of fundamental process 
simulation examples. In particular, the efficiency of the Pryce method is 
illustrated as a function of the number of tanks in process design. 

Keywords: differential algebraic equations, structural analysis, sym- 
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1 Introduction 

Differential-algebraic equations (DAE) systems arise naturally when modelling 
many dynamic systems. Dynamic process models and their properties form the 
background of any process control activity including model analysis, model pa- 
rameter and structure estimation, diagnosis, regulation or optimal control. In 
particular, the structural analysis of dynamic lumped process models forms an 
important step in the model building procedure pQ, and it is used for the deter- 
mination of the solvability properties of the model. Furthermore, the dynamic 
lumped process models often require the consistent initial conditions and solu- 
tion of high-index differential- algebraic systems. 

The index is a notion used in the theory of DAEs for measuring the distance 
from a DAE to its related ODE. High-index DAE systems need prolongation (dif- 
ferentiation) to reveal all the system's constraints, and to determine consistent 
initial conditions. The key steps include identifying all hidden constraints on for- 
mal power series solutions in the neighborhood of a given point, and are required 
to prepare the system for numerical integration. So for such differential systems, 
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prolongation is unavoidable. In the present work, Pryce developed a Taylor series 
method based on his structural analysis method (213) and on Pantelides' work 
in [4]. Pantelides' method gives a systematic way to reduce high-index systems 
of differential- algebraic equations to lower index, by selectively adding differen- 
tiated forms of the equations already present in the system. It is implemented 
in several significant equation-based simulation programs such as gPROMS [5], 
Modelica [B] and EMSO 7 . However, the algorithm can fail in some instances. 
Pryce's structural analysis is based on solving an assignment problem, which 
can be formulated as an integer linear programming problem. It finds all the 
constraints for a large class of ODE using only prolongation, which can be con- 
sidered as fast prolongation method. Corless et al. show Pryce's method can be 
extended to give a polynomial cost method for numerical solution of differential 
algebraic equations [5]. Wu et al. give a differential algebraic interpretation of 
Pryce's method for ODE, which generalizes to a certain class of PDE for finding 
missing constraints [5] . Mani shows how pre-symbolic simplification can usefully 
extend the applicability of the Pryce method on models produced by MaplcSim 

In |llll2j . Leitold et al. propose the structural analysis of process models 
using their representation graphs for the determination of the most important 
solvability property of lumped dynamic models: the differential index. Their 
graph-theoretical method depends on the change in the relative position of un- 
derspecified and overspecified subgraphs and has an effect to the value of the 
differential index for complex models. If these subgraphs move further from their 
original positions the value of differential index increases. In this paper, we con- 
sider other approaches for the structural analysis of dynamic lumped process 
models for high-index DAE systems. In particular, we consider Pryce's method 
and the symbolic differential elimination package rifsimp. Pryce's method is a 
robust and reliable method for remedying the drawback of the approach [11112] 
and doing so automatically. This is a powerful way to determine the index of the 
system, its number of degrees of freedom, and exactly which components should 
be given initial values. The key idea is taken from Pryce's signature-method. The 
nice feature of the work is a simple and straightforward method for analysing 
the structure of a differential algebraic system. 

The rest of this paper is organized as follows. Section 2 describes Pryce's 
method and introduces the symbolic differential elimination package rifsimp in 
Maple. Section 3 gives the structural analysis of simple process models using 
these approaches. Section 4 gives some experimental results. The final section 
concludes this paper. 



2 Preliminaries 



In this section, we give a brief review of Pryce's method and some remarks, and 
present the symbolic differential elimination with Maple's rifsimp package. 
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2.1 Pryce's method 



We review below the main steps of Pryce's structural analysis and the corre- 
sponding algorithm following [213] . We consider an input system of n equations 
/ = 0, where / = (/i, f% ■ ■ ■ , /„) in n dependent variables x\{t) , X2{t) , ■ ■ ■ ,x n (t). 
Step 1. Form the n x n signature matrix S = Oij of the DAE, where 

{highest order of derivative to which the variable 
Xj appears in equation fa, 
or — oo if the variable does not occur. 

Step 2. Solve an assignment problem to find a HVT (highest value transversal), 
which is a subset of indices describing just one element in each row and 

each column, such that ^ Uij is maximized and finite. 

Step 3. Determine the offsets of the problem, which are the vectors c = 
{ci)i<i< n ,d — (dj)i<j< n , the smallest such that dj — Cj > ov,, for all 1 < i < 
n, 1 < j < n with equality on the HVT. This problem can be formulated as an 
integer linear programming problem (LPP) in the variables c = (ci,C2, • • ■ ,c„) 
and d = (di,d 2 , ■ ■ ■ ,d n ): 



Minimize z = dj — a, (la) 

j i 

where dj — c, > <xy for all (lb) 

Q > for all i. (lc) 



The structural index is then defined as 

f for all dj > 
v = maxj a + < 

I 1 for some dj = 0. 

The structural index is no less that the differential index on first order DAE. 
Step 4. Form the n x n system Jacobian matrix J where 

j. . = f mjd^Jth derivative of Xj) ^ this derivative is present in ft 
1 otherwise. 

Step 5. Choose a consistent point. If J is non-singular at that point, then the 
solution can be computed with Taylor series or numerical homotopy continuation 
techniques in a neighborhood of that point. 



Remark 1. The computation of c and d only involves the information on dif- 
ferential order and is consequently very fast in Step 3 of Pryce's method. This 
problem is dual to the assignment problem. The time complexity of the assign- 
ment problem can be done at polynomial cost by using the Hungarian Method 
[T5] . Generally, such problems can be solved very efficiently in practice. 
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Remark 2. After we obtain the number of prolongation steps Ci for each equation 
from Step 3 of Pryce's method, we can enlarge the system of equations using 
c. We assume c\ > c-i > ■ ■■ > c„, and let k c = max^ Ci = c\, which is closely 
related to the index of DAEs. Consider the equations obtained by taking the 
^-derivative of = fx = up to the c^th derivative, 1 < i < n, that is the 
collection 

(ci) ^ 



AO) f (i) 

Jl iJl ; 



f(0) 



(Cn) 



= 0, 



(2) 



where W denotes d l /dt l 0. By the definition of cry and inequalities (lb), the 
derivatives of the Xj that occur in equations (2) all lie in this set: 



r (°) T « . . . T ( dl 



T (o) T (i) 



■ • ,X 



(d„) 



(3) 



Represent (3) as a vector X, then (2) can be written as a system 



= F(t,X) = 



F (t,X ) 
Fi{t, Xq, Xi) 



\ 



\Fk d (t,X ,Xi, - ■ ■ ,X kd -i,X kd )/ 



(4) 



where kd = maxj dj = di, and assume d% > d% > • • • > d n . In particular, 
for < i < k c , Fi has fewer variables than J^+i. The block structure form 
Bi(0 < i < k c ) in the case Cj = Cj+i + 1 is given in Table 1. 



Table 1: The triangular block structure of F for the case c, = c, + i + 1 



B 


B 1 




Bk c -i 


Bk c 




F w 

4 o) 


F (a) 


F (oi-i) 


F (cx) 

F (c 2 ) 



Remark 3. Fast prolongation produces a simplified system to which a standard 
numeric solver can be efficiently applied. 

1 W is defined by the same way for the rest of this paper. 
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2.2 Symbolic differential elimination 

Maple's rifsimp package can be used to simplify small- and middle-scale DAEs, 
and overdetermined systems of polynomially nonlinear PDEs or ODEs and in- 
equations to a more useful form [14]. For the DAEs and ODEs the only inde- 
pendent variable is time. It processes systems of polynomially nonlinear PDEs 
with dependent variables ttx, tt2, • • ■ ,u n , which can be functions of several inde- 
pendent variables. 

The key idea of algorithm is substitution and differential elimination, which 
requires a ranking to be defined on the dependent variables and their derivatives. 
A basic step of differential elimination algorithms linearly appearing is to write 
the system in solved form with respect to each highest ranked derivative. It 
is treated by methods involving a combination of Grober bases and Triangular 
decompositions. Another key step in such algorithms is the taking of integrability 
conditions between equations. 

The rifsimp algorithm is essentially an extension of the Gaussian elimination 
to DAEs and systems of nonlinear PDEs. It differentiates the leading nonlinear 
equations and then reduce them with respect to the leading linear equations. 
If zero is obtained, it means that the equation is a consequence of the leading 
linear equations. If not, it means that this equation is a new constraint to the 
system. This is repeated until no new constraints are found. See Section 3.3 for 
a simple example. 

3 Structural analysis of simple process models 

In this section, we apply the above techniques to structure analysis of dynamic 
lumped process models DAE systems. The model is taken from dynamic process 
simulation and multi-domain modeling and simulation of complex systems. Here, 
the cascade of perfectly stirred tank reactors yields the basic examples of the 
paper, see Fig.l. 

3.1 Main algorithm 

Suppose a system consists of k perfectly stirred tank reactor. A feed of con- 
centration Co(t) is fed into the first tank. The concentrations in the tanks are 
described by the following equation: 

Cf = ^(CViW-aW) 1 = 1,2,...,* (5) 

where Ci(t) is the concentration in the tank i, q(t) is the flow rates from tank to 
tank and Vi(t) is the fixed volume of the tank i. Thus the flow rates between the 
tanks qi(t) are all the same that qi{t) = q(t) — Q(t), where Q(t) is a specified 
function of t. 

In general, there are two different specifications that can be added to these 
equations according to the modelling goal: 

a) in dynamic simulation studies the feed concentration Co(t) is given by Co = 
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Fig. 1: Sequence of liquid tanks, where for the i-th tank and Fi(i = 

1,2, ••• , fc) are the inlet and outlet flow rate, Ci is the concentration and V% 
is the fixed volume of the i-th tank. 



C0(t); 

b) in dynamic design the product concentration Ck(t) is given by Ck(t) = Ck{t). 

When applied to process system a) and b), the main steps of our approach 
are: 

Step 1. Construct the original system as follows based on the equation §5§: 
F:= [d(t)W =q(t)(C i - 1 (t)-C i (t))/Vi(t),i = 1,2,- ■■ ,k, 

(6) 

V5(t)W =0, i = l,2,..., A, g(t)=Q(t)], 

where fc is the number of tanks. 

Step 2. Obtain the original condition and add it into F. There are two cases: 

a) in dynamic simulation the tank feed concentration Co(t) is given as a function 
of time then get 2fc + 1 equations in 2fc + 1 unknowns: Co(t) = C0(t), this 
is essentially index 1 no matter what fc is, and is a trivial system. Symbolic 
differential elimination can be used for case a); 

b) in dynamic design the product concentration Ck(t) is given as a function: 
Cfe(t) = Ck(t). It is a nontrivial system, which is high-index as fc increased. 

Step 3. Call the Pryce's algorithm of Section 2 to solve the vector c and 
d, and enlarge the initial system by fast prolongation. Alternatively symbolic 
differential elimination can be used for case b). 

Step 4. Check the Jacobian matrix J with the coefficients of highest deriva- 
tives equations and compute the consistent point. 

Remark 4- Based on the structure analysis of Pryce's method, it is practical and 
efficient for dynamic lumped process models DAE systems. In general, the goal 
of structural analysis of DAEs is to differentiate the equations in such a way 
that the coefficient (Jacobian) matrix of the highest derivatives is non-singular. 
It means that some equations need prolongations on independent variable to 
balance the coefficients matrix. So it can computing Jacobian matrix of the 
lower derivatives equations by an iterative procedure for finding all consistent 
points. 
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Remark 5. Compared with the structural analysis of process models using their 
representation graphs method, the advantages of our algorithm are: 

• We efficiently apply the fast numerical and symbolic computations to a wide 
variety of physical models generated by the equation-based technique. 

• For the large models, we can keep the structural index of system remaining 
unchanged. Moreover, the prolongation system has a favorable block triangular 
structure to compute the missing initial conditions more efficiently. 

3.2 Main results 

For the general dynamic lumped process models DAE systems, we can obtain 
the offsets of vector 



and 



c= (0,1, 2,.-- , fe-1, 0,0,1, 2,--. ,fc) 



(0,1,2,--- ,fc,l,l,2,--- , fc-1, fc-1) 



by Pryce's method. Therefore, we have the following ranking of dependent vari- 
ables. 



V k -i 

v k 
V 1 J 



-< 



(i) 



v k {1) 
V J 



-< 



c. 



(1) 



a 



(fc-i) 
fc 

Y-2 



r(fc-3) 
fc-1 



(fc-2) 
fc 



-< 



\ q (k-2) ) 



(2) 



a 



C (fc) 
(i) 



(1) 



(fc-2) 



fc-1 

(fc-1) 



\ q (k-i) J 



Based on the above ranking, we can obtain the sequence of solving initial 
value problem for the dynamic lumped process models DAE systems. It is equiv- 
alently the block-triangular system that has full row rank for each fc. 

From ©, F and the original condition C k {t) = Ck(t) have 



fc-i 



fc 



M = (J2 Ci) + (2k + 2) = (Q2 i) + {J2 0) + (2fc + 2) = fc 2 + 2fc + 2 



l l 

components. The number of variables is 

fc fc-i 



N = (5^) + ( 2fc + 2 ) = ((5Z^ + 1 + E- ? ') + ( fc_1 )) + ( 2fc + 2 ) = fc2 + 3fc + 2. 
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Considered as M algebraic equations in N variables, it has a solution (t*,X*), 
where 

and that J is non-singular. Therefore, (t*,X*) is a consistent point. In a neigh- 
borhood of the point (t*,X*), the solution manifold has D degrees of freedom 
®. 

Lemma 1. At a point (t*,X*) in M where J is non-singular, M is locally a 
manifold of dimension tt + 1 parameterized. The solution manifold has D degrees 
of freedom, where 

D = tt = ^2dj - 53 Ci = N - M. 

The above shows that if we find a solution (t*,X*), this is a consistent point, 
and if the number of degrees of freedom D > there are other consistent points 
nearby for the same t. 

Theorem 1. The general dynamic lumped process models DAE systems have 
degrees of freedom D = ^ dj — ^ Cj = k, where k is the number of tanks. The 
structural index is k + 1 . 

Proof. From Lemma 1 , we have degrees of freedom 

k k-i 



i i 

fc-l k-2 

((E 3) + (E j) + k - 1 + k ) = k + 1 + k - 1 + ° - k = k - 



(7) 



Because the d\ = 0, the structural index is max^ (ci) + 1 = k + 1. 

Here, we give the degrees of freedom and structural index of the general dynamic 
lumped process models DAE systems that is a function of the number of tanks 
k. 



3.3 A detailed example 

Example 1. We propose a simple example to set the number of tanks k := 3 case 
a) to illustrate the rifsimp algorithm. 

Step 1: Construct the original system as follows: 

(1) g (t)*(C (t)-Ci(t)) (1) g(t) * (C^t) - C 2 (t)) 
sys := [<*(*)< = ^ ,C 2 (tY = ^ , 

Cs {t)W = q(t)*(C2{t) -C 3 (t)) t ^ (t)(1) = Qj Va(t)(1) = Qj m = Q) ff(t) = g(t)] 

V3W 
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Step 2: Obtain the original condition Co(t) = C0(t), and add it to sys; 
Step 3: Call rifsimp algorithm to reduce the system as follows: 

[Ci(t) = W) ' 2( } = W) ' 

(1) = g W ^ 2 (t) 7 Qft),c 3W ^ i(t)a) ^ 0it/2(t)(1) ^ 0; 

V3{t) 

t/ 3 (i) (1) = 0,C o (t) = C0(t),g(t) = Q(t), ^ 0, V 2 (t) ^ 0,^ 3 (*) ^ 0]. 

Remark 6. In this paper, we consider the modelling goal for case a) by the 
rifsimp algorithm. The main reason is the specific structure of models, which 
is the quasi-triangular system and has Co(t) = C0(t) specified. Therefore, it is 
only simple check. But it becomes rapidly more complicated as the number k 
increased for case b). 

Example 2. We propose a simple example to set the number of tanks k := 4 case 
b) and illustrate our algorithms. 

Step 1: Construct the original system as follows: 

sys := [Di = C,( t )<« - '"»• taW -C-W) = , D2 = a(() (« 

w = °- 3 = ) W) = ' 4 = 4 ) 

_ g (t)*(c,(t)-c,(t)) = Q; Ds = (1) = Q> D6 = (1) = Qj 

#7 = V 3 (t)W =0,D 8 = V 4 (t)W = 0,D 9 = q(t) - Q(t) = 0]; 

Step 2: Obtain the original condition Ci{t) = C4(t), and add D w = 

Ci(t) - C4(i) = to sys; 
Step 3: Obtain the variables list variables := [Co, Ci, C2, C3, C4, V\, V2, V3, V4, g]; 
Step 4: Call the Pryce's method and solving this integer LPP by LPSolve 

in the Optimization package of Maple, we obtain the fast prolongation 

times for the i-th equation from c, and the highest order of derivative 

variables from d as follows: 

ci = 0, c 2 = 1, c 3 = 2, c 4 = 3, c 5 = 0, c 6 = 0, c 7 = 1, c 8 = 2, c 9 = 3, c i0 = 
4, 

g?i = 0, d 2 = l,e?3 = 2,^4 = 3,c? 5 = 4,c? 6 = l,c?7 = l,dg = 2,c? g = 
3, dio = 3. 

Therefore, according to the Cj values it can be prolonged for the corre- 
sponding equations automatically. Enlarged sets of variables: 

i^Oi <^1> W 1^2,^2 i L/ 2 i°3,^3 )<^3 ) 3 i L/ 4,L'4 ,^4 ,^4 ,^4 , 

equations: 

{D 1 ; D 2 , ; D 3 , , ; D 4 , , , D< 3) ;D 5 ;D 6 ;D 7 , ; D 8 , £>W , ; 



10 



Xiaolin Qin et al. 



D 9 ,D£\D { g\D^ ) ;D w ,D§,D<$,D[ ! $,D[$}. The system Jacobian J 
is: 



J := 



Vi(t) 



q(*) 

V2(t) 











v 3 (t) 

v 4 (*) 
o 

























q(t) (djt) - C 2 (t)) 







1 

10 














q(t) (C 2 (t) - C 3 (t)) o 



q(t) (C 3 (t) - C 4 (t)) 



C 3 (t)-Ci(t) 



Vk(t) 




1 






Step 5: Computing the Jacobian matrix J = — ^ ^ ^ ^v 3 ^y 4 ( t y which 
is non-singular. And then we can compute the consistent point by nu- 
merical methods, such as Taylor series methods, Homotopy methods. 



Remark 7. In particular, we can obtain the coefficient (Jacobian) matrix that 
is sparse dramatically. For the highest derivatives, the determinant of Jacobian 
matrix is det J — — Vi(t)V2(t)—v k (t) wnere & i s the number of tanks. 



4 Experimental Results 



An efficient practical implementation of Pryce's method is in Maple. The follow- 
ing examples run in the platform of Maple and Intcr(R) Corc(TM) i3 2.40GHz, 
2.00G RAM. We give some experimental results using symbolic differential elim- 
ination and fast prolongation for structural analysis of dynamic lumped process 
models DAE systems. In Fig. 2, we present the time for symbolic differential 
elimination by rifsimp of Maple package and fast prolongation as the number of 
tank reactors k increased. In Fig. 3, we present the memory usage for symbolic 
differential elimination by rifsimp of Maple package and fast prolongation as the 
number of tank reactors k increased. 

The system Jacobian is very sparse for case b). Its determinants are evaluated 
symbolically to be det J = — y 1 ( t )v|(t)...v fc (t) where k is the number of tanks, and 
is non zero. In other examples the alternative is to usually find an approximate 
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point satisfying the constraints by numerical method (eg. Homotopy method) 
and evaluate the condition number of the Jacobian to carry out he validation. 




2 3 4 5 | 7 8 9 10 11 12 

k (number of tanks) 



Fig. 2: Time for structural analysis of dynamic lumped process models DAE 
systems using symbolic differential elimination and fast prolongation. 




k (number of tanks) 



Fig. 3: Memory usage for structural analysis of dynamic lumped process models 
DAE systems using symbolic differential elimination and fast prolongation. 



From Figures 2 and 3: 
• The time of structural analysis of dynamic lumped process models DAE sys- 
tems for fast prolongation is small and ultimately grows slowly in the range 
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of degrees of freedom considered. The time for symbolic differential elimination 
method grows much faster. The main reason is that fast prolongation only needs 
to solve an integer linear programming problem, but the symbolic differential 
elimination needs a large number of eliminations and differentiates. Therefore, 
the symbolic differential elimination is more difficult for the general high-index 
DAE systems. 

• The memory shows steady growth as the number k increases. The memory 
usage of symbolic differential elimination grows very quickly. 

The above analysis and experimental results, motivates consideration of hy- 
brid techniques involving a combination of symbolic differential elimination and 
fast prolongation for large DAE models. However, symbolic computations have 
the disadvantage of intermediate expression swell. In the future, we would like to 
consider a combination of partial symbolic differential elimination and fast pro- 
longation to model and simulate realistic physical models. We hope to give the 
specific structural analysis algorithms that exploit the form of systems appearing 
in applications. 

5 Conclusion 

In this paper, we have investigated the high-index structural analysis problem 
for the class of dynamic lumped process models DAE systems by Pryce's method 
and symbolic differential elimination. We designed the algorithm to automati- 
cally analysis the structural of simple process models, and showed that the rif- 
simp algorithm of Maple package reduces the original system to standard form. 
We also gave the degrees of freedom and structural index of the dynamic lumped 
process models DAE systems that is a function of the number of tanks k. More- 
over, those approached can be generalized to a wide variety of physical models 
and analyzed the structural of square and non-square nonlinear DAE and PDAE 
systems. 
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